# Replication archive for: Alexander Coppock, Kimberly Gross, Ethan Porter, Emily Thorson, and Thomas J. Wood, 
# "Conceptual Replication of Four Key Findings about Factual Corrections and Misinformation during the 2020 US Election: Evidence from Panel-Survey Experiments." 
# Forthcoming in the British Journal of Political Science doi:10.1017/S0007123422000631

library(tidyverse)
library(coefplot)
library(estimatr)
library(nnet)
library(lmtest)

all_panels_long <- read_rds("data/all_panels_long.rds")


all_panels_long <-
  all_panels_long %>%
  mutate(
    race_asian = as.numeric(lucid_race == "Asian"),
    race_black = as.numeric(lucid_race == "Black"),
    race_hispanic = as.numeric(lucid_race == "Hispanic"),
    race_other = as.numeric(lucid_race == "Other"),
    race_white = as.numeric(lucid_race == "White"),
  )



balance_regs <-
  all_panels_long %>%
  group_by(fc, panel_factor) %>%
  do(tidy(lm_robust(
    cbind(
      lucid_pid_7n,
      lucid_age,
      race_asian,
      race_black,
      race_hispanic,
      race_other,
      race_white,
      lucid_hhin,
      political_knowledge_pre,
      political_interest_pre,
      cognitive_reflection_pre,
      need_for_cognition_pre,
      extraversion_pre,
      agreeableness_pre,
      conscientiousness_pre,
      emotional_stability_pre,
      openness_to_experience_pre
    )
    ~ treatment,
    data = .
  )))



balance_regs <-
  balance_regs %>%
  filter(term != "(Intercept)") %>%
  mutate(
    term = str_remove(term, "treatment"),
    cov_set = case_when(
      outcome %in% c(
        "race_asian",
        "race_black",
        "race_hispanic",
        "race_other",
        "race_white"
      ) ~ 1,
      outcome %in% c(
        "lucid_pid_7n",
        "lucid_age",
        "lucid_hhin",
        "political_knowledge_pre",
        "political_interest_pre",
        "cognitive_reflection_pre",
        "need_for_cognition_pre"
      ) ~ 2,
      outcome %in% c(
        "extraversion_pre",
        "agreeableness_pre",
        "conscientiousness_pre",
        "emotional_stability_pre",
        "openness_to_experience_pre"
      ) ~ 3
      
    )
  )


pos = position_dodgev(height = 0.5)

g <-
  ggplot(data = NULL) +
  aes(
    estimate,
    y = term,
    group = fc,
    shape = (p.value <= 0.05),
    color = (p.value <= 0.05)
  ) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_point(position = pos) +
  geom_linerange(aes(xmin = conf.low, xmax = conf.high),
                 position = pos) +
  scale_color_manual(values = c(gray(0.5), gray(0.1))) +
  facet_grid(panel_factor ~ outcome, scales = "free_x") +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.text.x = element_text(size = 6),
        axis.title.y = element_blank(),
        legend.position = "none") +
  labs(x = "Estimated effect on pre-treatment covariate, relative to the misinformation condition")

g1 <- g %+% filter(balance_regs, cov_set == 1)
g2 <- g %+% filter(balance_regs, cov_set == 2)
g3 <- g %+% filter(balance_regs, cov_set == 3)

# ggsave(filename = "plots/balance1.pdf", 
#        g1, width = 8, height = 9)
# ggsave(filename = "plots/balance2.pdf", 
#        g2, width = 8, height = 9)
# ggsave(filename = "plots/balance3.pdf", 
#        g3, width = 8, height = 9)


# In text calculations
sum(balance_regs$p.value <= 0.05)
sum(p.adjust(balance_regs$p.value, method = "BH") <= 0.05)
42 / 816


# Joint Test ------------------------------------------------

data_splits <-
  all_panels_long %>%
  split(list(.$fc, .$panel_factor))
     
joint_test <- rep(NA, length(data_splits))

for(i in 1:length(data_splits)){
  dat <- data_splits[[i]]
  fit <- multinom(
    treatment ~ lucid_pid_7n +
      lucid_age +
      lucid_race +
      lucid_hhin +
      political_knowledge_pre +
      political_interest_pre +
      cognitive_reflection_pre +
      need_for_cognition_pre +
      extraversion_pre +
      agreeableness_pre +
      conscientiousness_pre +
      emotional_stability_pre +
      openness_to_experience_pre,
    data = dat
  )
  joint_test[i] <- lrtest(fit)$`Pr(>Chisq)`[2]
}

gg_df <-
  tibble(raw = joint_test) %>%
  mutate(adjusted = p.adjust(raw, method = "BH")) %>%
  pivot_longer(cols = c(raw, adjusted)) %>% 
  mutate(name = factor(name, levels = c("raw", "adjusted")))



g <-
  ggplot(gg_df) +
  aes(value) +
  geom_histogram(binwidth = 0.05) +
  facet_wrap(~name) +
  theme_bw() +
  theme(strip.background = element_blank()) +
  labs(x = "p.value of joint test")
g
# ggsave("plots/joint_balance.pdf", g, width = 6.5, height = 4)


